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The generation of higher harmonics of the magnetoplasmon frequency which has recently been 
reported in strongly coupled two-dimensional Yukawa systems is investigated in detail and, in ad- 
dition, extended to two-dimensional Coulomb systems. We observe higher harmonics over a much 
larger frequency range than before and compare the theoretical prediction with the simulations. 
The influence of the coupling, structure, and thermal energy on the excitation of these modes is 
examined in detail. We also report on the effect of friction on the mode spectra to make predictions 
about the experimental observability of this new effect. 
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I. INTRODUCTION 

The behaviour of two-dimensional (2D) systems is of 
continuing interest in many fields of physics as the re- 
duced number of dimensions give rise to a number of pe- 
culiar properties. In highly correlated systems, where the 
potential energy due to the interaction dominates over 
the thermal energy (for recent overviews see [1, 2 ), 2D 
many-particle systems exhibit a strongly collective be- 
haviour which manifests itself, e.g., in anomalous trans- 
port properties of 2D liquids [3Hl]. Correlational effects 
are also responsible for additional shear mode excitations 
of 2D liquids O [9] which do not occur in weakly cou- 
pled systems but have been experimentally observed in 
strongly coupled dusty plasmas [TOl [11] . 

Subjecting a 2D many-particle system to a perpendicu- 
lar magnetic field gives rise to yet another line of research 
into the physics of low-dimensional systems [12^. The 
magnetic field effectively "mixes" longitudinal and trans- 
verse excitations, leading to two hybrid modes, the mag- 
netoplasmon and the magnetoshear, which have in recent 
times been studied in Coulomb and Yukawa systems [13]- 
[15] . These modes are well understood from a theoretical 
perspective, including descriptions based on the quasi- 
localized charge approximation (QLCA) [161 [IZ] oi" 
harmonic lattice approximations [9i il5j . 

However, besides these established modes, 2D Yukawa 
systems at strong coupling have recently been found to 
support additional high-frequency modes. These modes 
appear as higher harmonics of the magnetoplasmon in 
non-dissipative 2D Yukawa systems [18] and are reminis- 
cent of the classical Bernstein modes [T9]. Unlike these, 
however, the observed higher harmonics are not a pure 
magnetic effect but are, additionally, fundamentally af- 
fected by the strong correlations between particles, as 
was demonstrated in Ref. [TSl . 

In this work, we put the theoretical predictions of our 
earlier work [18 to the test by using more detailed simu- 



lations and observing considerably more high-frequency 
modes. In addition, we include the Coulomb case of van- 
ishing screening in our analysis, and demonstrate that 
the higher harmonics are also generated under these cir- 
cumstances, which are not only quantitatively but also 
qualitatively different from screened-interaction systems. 
Finally, we assess the possibility of experimental verifica- 
tion by investigating the relative intensities of the modes 
and including dissipative effects in our simulations. 

The remainder of this article is structured as follows: 
In Section [ll| we introduce our model and give details of 
the numerical procedure and the system of units. Sec- 



tion III introduces the longitudinal and transverse fluc- 
tuation spectra and the means by which they are com- 
puted. The results for non-dissipative and dissipative 



systems are presented in Section IV before we summa- 
rize our findings in Section [Vj 



II. MODEL AND DETAILS OF THE 
NUMERICAL SIMULATION 

Our model system consists of N point-particles situ- 
ated in a quadratic simulation box in the x, ?/-plane with 
side- length L. The particles are uniform in mass m and 
charge q and are subject to periodic boundary conditions 
to avoid surface effects. The particles propagate accord- 
ing to the coupled equations of motion. 



mvi = Fi ^ qr-i X B ^ Si, 



(1) 



where the force Fi follows from the Yukawa potential 
created by all other particles. 
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Here, Xd denotes the Debye screening length and the 
primed sum indicates the omission of the term j = i. 
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For Xd oo, we recover the weh-known one-component 
plasma (OCP). The magnetic field is oriented perpen- 
dicular to the plane of the particles, B = Be^^ and the 
Langevin term Si in Eq. ([T]) is defined as 

Si = -rriivri ^ Ri. (3) 

This Langevin term is only included in the simulations 
with friction and Riit) is a Gaussian white noise with 
zero mean and the standard deviation 

{RaAto)R^,j{to + t)) = 2kBTiy6,j5a^ 5{t), (4) 

where a, /3 G {x, y} and T is the temperature. 

Eq. ([T]) is solved simultaneously for N = 4080 particles 
using standard molecular dynamics simulation with inte- 
grators adopted to the influence of the magnetic field [20j . 
When friction is included in the simulations, we apply 
an additional Ornstein-Uhlenbeck process in momentum 
space [m [22]. For Coulomb systems, the appropriate 
Ewald summation techniques to calculate the forces are 
employed [23l[24]. Starting from a random configuration 
of particles, the system is brought into equilibrium by a 
repeated rescaling of the particles' momenta according to 
the desired temperature T. After equilibrium conditions 
are realized, the system is advanced only according to 
Eq. 

In the following, lengths are given in units of the 

]^ /2 

Wigner-Seitz radius a = [nir] ' , where n is the areal 
number density of the particles. Time is given in mul- 
tiples of the inverse of the nominal angular plasma fre- 

1/2 

quency cjp = [ng^/ (2^0 ma)] . The system is charac- 
terized by four dimensionless parameters: 

1. the temperature T of the particles, given in terms 
of the system's nominal (i.e.. Coulomb-) coupling 
parameter F = / {A-KSQakBT)^ 

2. the inverse of the Debye screening length Xd-, nor- 
malized by the Wigner-Seitz radius, n = a/Xjj^ 

3. the strength of the magnetic field, expressed as the 
ratio P between the angular cyclotron frequency 
ujc = qB/m and the nominal angular plasma fre- 
quency, p = UJc/uJp, 

4. the strength of friction, given by v = u/ujp. 

III. COLLECTIVE MODES 

The collective excitations of the strongly coupled, in- 
teracting particles can be analysed through the micro- 
scopic excitation spectra. The dispersion, i.e., the re- 
lation between the wavelength and the wave frequency 
of the collective excitations, can be obtained from the 
analysis of the autocorrelation (ACF) of the density and 
current fiuctuations. For the Fourier components of the 
density fiuctuations, we have [25j 

N 

p{k,t) = ^e^'^-^-^*). (5) 



The dynamical structure factor S{k^uj) follows as (as- 
suming k = kcx and dropping the vector notation for k): 

S{k,uj) = ^^lim (6) 

where J^t denotes the temporal Fourier transform. 
The current operator is given by j25l [26] 

N 

j{k,t) = ^Vj{t)exp[ikxj{t)]. (7) 

i=i 

Separating longitudinal [X{k,t)] and transverse [r(/c,t)] 
currents, one arrives at the microscopic quantities [27] 

N 

Vjxex.p[ikxj], (8) 

N 

r{k^t) — Vjy exp[z/cxj], (9) 

from which one can then obtain the fiuctuation spec- 
tra L{k^uj) and T(k^uj) analogously to Eq. ([6| (replacing 
p{k^t) by X{k^t) and r{k^t), respectively): 

T(fc,w) = lim ^\Mrik,t)}\\ (11) 

The dynamic structure factor S{k^uj) can be shown to 
be related to L{k^uo) as [28] 

S{k,uj) = ^L{k,uj). (12) 

S{k^uj)^ therefore, contains no additional information 
over L{k^uo)^ so we concentrate on L(k^uj) and T{k^uj) 
in the following [29 . 

Our main interest lies in the mode spectra in the high- 
frequency range. We numerically evaluate the micro- 
scopic fiuctuations ([8| and ([9| at multiples of the min- 
imum wavevector k^i^ = {27r/L)ex as dictated by our 
use of periodic boundary conditions. The data is subse- 
quently Fourier analysed and the periodogram estimate 
of the power spectrum density is computed to obtain 
L{k^uj) and T{k,uj). Our simulations are typically ther- 
malized for 20.000 plasma periods and data are collected 
during 400.000 or more plasma periods. The spectra are 
computed for different values and combinations of the 
coupling r, screening magnetic field f3 and friction v. 

IV. RESULTS 

In this section, we study the generation of the higher 
harmonics in detail. Using more accurate data, we con- 
firm the theoretical formula of the mode spacing reported 
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in Ref. [TF in section IV A and then discuss the effect of 
the couphng, the structural properties and the thermal 
energy on the higher harmonics in frictionless systems in 
section [IVB| In section [TV C| we include dissipation and 
stochastic noise in our simulations to gauge the prospect 
of experimental verification. 



A. Dependence of the spectrum on /3 

To obtain a first overview of the collective excitations, 
we depict in Fig. [l] several results for the longitudinal and 
transverse wave spectra at different strengths of the mag- 
netic field for a Yukawa system with k = 2 and F = 200. 
The density plots in Fig. [l] show the sum of transverse 
and longitudinal excitations. For the unmagnetized case 
(top row), the mode-spectrum is well-known and accu- 
rately described, e.g., within the QLCA [8^. The two 
branches associated with strongly coupled, unmagnetized 
Yukawa and Coulomb liquids (transverse and longitudi- 
nal branch) are not discernible on the scales of Fig. [l] At 
higher frequencies, the mode spectra appear structureless 
(cf. right columns of Fig. [T]). 

In the the magnetized case (lower three rows of Fig. [l]), 
two features are immediately noticeable: i) The two 
modes from the unmagnetized case are now replaced by 
the upper- and lower-hybrid modes (UH and LH) with 



the ^ limits uj^^ = and uj^^ = P^- h) The quali- 
tatively new feature reported in Ref. [18] is the emergence 
of multiple new branches at higher frequencies. These 
are higher harmonics reminiscent of the Bernstein modes 
which have been observed half a century ago in strongly 
magnetized ideal plasmas [19]. Typical for these modes 
is that their frequencies are equally spaced, appearing at 
multiples of the cyclotron frequencies ujc and that they 
are undamped and exist in the entire wave number range. 

The present modes show a number of fundamental 
differences which are caused by correlation effects: all 
modes are damped, and they appear only beyond a criti- 
cal wave number which increases with the mode number. 
The most important difference is the mode frequency. In- 
stead of being spaced by the cyclotron frequency cJc, here 
the modes appear at multiples of the magnetoplasmon 
frequency. The frequency of these "dressed Bernstein" 
modes, cj^, n = 2, 3, . . . , were found to be well described 
by the relation [18 



l,oo' 



Lol^iP,K)=u;^,iP) + 2u;%iK). (13) 



Here, the dominant single-particle oscillation frequency 
is denoted by uje^ the Einstein frequency, which strongly 
varies with hz but also weakly depends on F [30 . Some 
values of a;£;(r, /^) are collected in Tab. [l[ 

In Ref. [18 , the data allowed to detect higher harmon- 
ics of the magnetoplasmon up to the third order [(n = 4 




FIG. 2: (Color) Combined mode spectrum L{k^uj) -\-T{k,(jj) 
for a) At = 0.0, r 100 and h) k = 2.0, F 200, and /3 1.0. 
Black lines are the theoretical prediction for higher harmonics 
according to Eq. (13). 



in Eq. (13)]. The present data show higher harmonics up 
to n = Tjcf. Fig.jl] third row). This provides us with the 
opportunity to verify the validity of Eq. ( 13 ) to a much 
higher precision. To this end, we depict in Fig. |2] the 
combined mode spectra for /3 = 1.0 in two different sys- 
tems. The solid lines indicate the theoretical prediction 



of Eq. (13). A very good agreement between the the- 
ory and the simulations is evident. In addition, Fig. |2|a) 
illustrates that the generation of the dressed Bernstein 
modes also occurs in a Coulomb system in which the 
particle interaction is long-ranged. This feature substan- 
tially expands the scope of the higher harmonics gener- 
ation to other fields in plasma physics, including, e.g., 
ionic plasmas in traps [31^ and electron-hole plasmas in 
semiconductors. 

We now quantify the higher harmonics in more detail. 
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FIG. 3: K = 2.0, r = 200, and k = 0.0, F = 100 at /3 = 1.0, 
ka = 2.0 a) The peak position of the nth mode. The straight 
lines are the theoretical predictions |T3| with (jJe/(jOo — 0.62 
[hv — 0.0) and ixJe/(jOq — 0.34 [hi — 2.0). h) The peak intensity 
normalized to the intensity of the magnetoplasmon (n = 1). 
An exponential decay is shown by the straight line. 



Fig. [3] shows the frequency and oscillator strengths of the 
higher harmonics and the magnetoplasmon at ka = 2.0 



and fixed F and (3. The agreement between relation (13) 



and the peak position from the simulations is again ex- 
cellent (Fig.[3|i), both for the OCP and a typical Yukawa 
system with = 2. Note that no free parameter enters 
the theoretical prediction and that the differences in the 
theoretical predictions are only mediated by uje- 



K 


F 







100 


0.62 


1 


150 


0.52 


2 


40 


0.40 


2 


100 


0.36 


2 


200 


0.34 


3 


600 


0.20 



TABLE I: The Einstein frequencies ue d^s calculated from 
QLCA |2j for various combinations of F and 



A quantity of central interest is the relative intensity 
of the generated higher harmonics. These values provide 
one with a first estimate of the required experimental 
sensitivity to observe the described effects. Figure ^ 
shows the relative intensities at ka = 2.0 for the Yukawa 
and Coulomb spectra. Evidently, the data points are 
well described by a decay of the form 10~2-^(^~^). The 
first of the higher harmonics (i.e., n = 2) appears with an 
intensity of about one hundredth of the magnetoplasmon. 
Such intensities are observable against the background 
noise in typical dusty plasma experiments [TOj . 
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FIG. 4: (Color) /3 = 0.7, n = 2.0, T = 10, 40, 200, 1000, (from top to bottom) Left: Collective excitation spectra, L{k, uj)+T(k, u), 
Right: L{k^ uj) + T(/c, uj) at two fixed values of ka — 2.0; 5.0 and values of V as before. 



B. Dependence of the spectrum on the interaction, 
temperature, and structural order 

The frictionless Yukawa system is (at fixed magneti- 
zation p) characterized by two parameters: The inverse 
temperature (Coulomb coupling parameter) T and the 
interaction range (inverse Debye length) n. The influ- 
ence of the temperature is investigated by changing, at a 
fixed value of the coupling parameter V. 

In Fig. [4j different situations are depicted for a sys- 
tem with K = 2^ ranging from intermediate coupling 
(F = 10; 40) to strong (F = 200) and very strong cou- 
pling (F = 1000). Note that for F = 1000, the sys- 
tem is already far in the microcrystalline regime, F^ei = 
r/Fmeiting ^ 2.4. The decreased thermal background 
makes the fundamental and higher harmonics stand out 
very clearly at strong coupling, and the signal-to-noise 
ratio is higher in these cases. For example, harmonics 
with n > 3 are far easier to distinguish at F = 1000 
than at F = 40. The intensity ratio M = In/Io between 
the intensity of the nth harmonic and the fundamental 
is, however, clearly decreased for stronger coupling [32 . 
This indicates that a certain degree of disorder is required 
for an effective generation of higher harmonics, which in 
turn might depend on anharmonic effects that are less 
important in highly ordered systems. 

Having established the influence of the temperature 



(inverse coupling), we now determine the influence of the 
interaction range. To this end, in Fig. [Sj we show the 
combined mode spectra at a fixed F = 100 for different 
values of hz and two values of ka. The curves appear very 
similar, with comparable signal-to-noise-ratios and sim- 
ilar relative intensities of the higher harmonics. This is 
despite the strong differences in structural order among 
the systems reflected in the different radial pair distribu- 
tion functions (RPDF, inset of Fig[5| defined as 



9{\r\) 




Y^S{\r-r,,\) 



(14) 



The sole effect of an increase in on the higher har- 
monics is a systematic shift toward smaller frequencies. 
This result is not unexpected and indeed predicted by 
Eq. (13): With increasing /^, the particle oscillations de- 



crease in frequency as the interparticle potential becomes 
more lenient. The Einstein frequency, therefore, dimin- 
ishes with increasing k, (cf. Table |l| and, consequently, so 
do the frequencies at which the higher harmonics appear. 

So far we have investigated the influence of F and hi on 
the higher harmonics. A change of either of these param- 
eters also affects the structural properties of the system. 
By changing F and hz simultaneously in such a way that 
the RPDF remains practically unchanged [30l f33H35] , we 
can compare systems with identical structural proper- 
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K = 1.0 Trei = 0.73 K = 3.0 r^ei = 0.08 




2 4 6 8 10 12 

uj/ujp 



FIG. 5: (Color) Combined mode spectra at ka = 2.0; 5.0 for 
identical temperature (F = 100) and k = 1.0; 2.0; 3.0 (highest 
to lowest curve), /3 = 1.0. The curves are shifted vertically 
for clarity. The inset shows the RPDF g{r). 

ties. We achieve this by fixing the ratio of the tempera- 
ture to the point of (micro-)crystallization, F/Fmeiting = 
Frei [36 , 37 . Systems corresponding to the same F^ei 
but different F have thus the same structural order but 
different thermal energy. 

The results are presented in Fig. [6j where a fixed value 
of Frei is maintained for systems with different values 
of K,. A combination of the two previously described 
effects is observed: The higher harmonics shift toward 
lower frequency with increasing tz^ and decay more slowly 
albeit with a decreased signal-to-noise-ratio with increas- 
ing temperature (inverse coupling). 

In conclusion, we find that the generation of higher 
harmonics is nearly independent of the intricate details of 
the structure and the range of the interaction between the 
particles (except for a systematic frequency shift) but, 
instead, is dominated by the thermal energy available to 
the system. 



C. Influence of dissipation on the higher 
harmonics. Prospects for experimental observation 

Besides noise resulting from the measurement itself, 
many physical systems are also influenced by fluctua- 
tions due to the coupling to their surroundings. In many 
cases, these fluctuations and dissipation effects are well 
described by a constant friction and stochastic noise, an 
approach which has been used extensively in the descrip- 
tion of, e.g., dusty plasmas. We now introduce dissipa- 
tion and stochastic noise in our MD simulations to esti- 
mate the effect of friction-induced noise on the fluctua- 



K = 1.0T= 90 /^ = 3.0r = 600 




uj/ujp 

FIG. 6: (Color) Combined mode spectra at ka = 2.0; 5.0 for 
identical structural properties corresponding to Frei = 0.48 
but different k, = 1.0; 2.0; 3.0 (highest to lowest curve), /3 — 
1.0. The inset shows the RPDF g{r). 

z/ = 10-1 z/ = 10-3 

z/ = 10-2 z/ = 10-^ 
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FIG. 7: F = 200, k = 2.0, /3 = 0.5: Combined mode spectra 
for different values of the friction strength, from u — 10~^ 
(top curve) to v — 10~^ (lowest curve). 



tion spectra. In other words, we now include the term ^ 
in Eq. 

In Fig. [7j we depict the longitudinal fluctuation spec- 
trum L{k^uj) for fixed values of /c, /3 and F, while varying 
the magnitude of the friction u. As expected, the inclu- 
sion of friction has a two- fold effect: a) The overall noise 
level in the spectrum is increased and begins to overlap 
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the higher harmonics and b) the peaks in the spectra 
become increasingly broader with higher friction. 

The noise resulting from the friction clearly dominates 
the intrinsic noise for the data depicted in Fig. [Tj The 
high-cj limit of the spectra increases linearly with the 
friction coefficient. From the data presented in Fig. [7| 
we conclude that at a moderately coupled liquid state, 
F = 200 and k = 2, the detection of the second harmonic 
is possible at /3 = 0.5 and a friction of u ^ 0.01 (note that 
u is given here in units of the nominal plasma frequency) . 

We now give estimates for the plasma parameters re- 
quired to observe the generation of higher harmonics un- 
der experimental conditions. The main obstacle in at- 
taining our simulation conditions in experiments are the 
required high magnetic fields. Combining the formulas 
for the cyclotron frequency uJc and the plasma frequency 
oOp and solving for the magnetic field, we obtain (in SI 
units): 

where p is the mass density. Typical 2D dusty plasma 
parameters are R = 1 . . . 10/im, p ^ Ig/cm^ and a = 
0.1... 1mm, giving rise (for /3 ^ 0.074) to a required 
magnetic field strength of about 20T, which is at the edge 
of current experimental possibilities. That number can 
be reduced by using lighter particles or compressing the 
system to decrease the inter-particle spacing (an increase 
of the average density by a factor of, e.g., 2 lowers the 
magnetic field requirements by 40%). Also note that the 
required magnetic field is independent of particle charge, 
which, however, enters into the coupling parameter F. 

The dissipation in 2D dusty plasma experiments can be 
quite low, a typical value is u ^ 0.02 [3 . Thus, according 
to our simulations, friction at this rate should not prevent 
the generation of the higher harmonics. It is, however, 
not uncommon to experience friction rates in excess of 

= 0.1, at which the observation of the higher harmonics 
is not possible. In conclusion, the observation of higher 
harmonics appears as a possible experimental venue but 
poses several challenging restraints on the experimental 
setup. 

V. SUMMARY 

In this work, we have compared the previously derived 
theoretical description of the mode spacing of the higher 



harmonics with extended MD simulations and found ex- 
cellent agreement over a wide range of coupling strengths 
and interaction range. We have also, for the first time, 
demonstrated that the generation of higher harmonics is 
present in strongly correlated systems with Coulomb in- 
teraction as well, which is especially interesting in the 
light of additional possibilities for experimental verifica- 
tion. The relative intensities of the higher harmonics 
have also been investigated and were found to decay, to 
a very good approximation, exponentially with the order 
of the harmonics. Finally, we have included dissipation 
in our simulations and estimated the maximum accept- 
able level of such friction for experiments. Combining 
estimates for the strength of the magnetic field and the 
maximum friction levels, we conclude that the experi- 
mental observation of higher harmonics is possible but 
challenging. 

Finally, we expect that our results are also of direct 
relevance for strongly correlated quantum plasmas. In 
the case of magnetized ideal quantum plasmas the "tra- 
ditional" Bernstein modes have been predicted theoret- 
ically long ago by Horing et al. |38] and were observed 
experimentally in electron- hole plasmas, e.g. [39^. In 
recent years strongly correlated quantum Coulomb sys- 
tems, including liquid states and Wigner crystals of elec- 
trons [40 , holes [41 and excitons, e.g. [42 , moved into 
the focus of research. It is expected that these systems, 
if placed into a strong magnetic field, should exhibit a 
similar spectrum of "dressed Bernstein modes" as their 
classical counterpart studied in the present paper. 
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